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ABSTRACT 

The merger of two white dwarfs (WDs) creates a differentially rotating remnant which 
is unstable to magnetohydrodynamic instabilities. These instabilities can lead to vis- 
cous evolution on a time-scale short compared to the thermal evolution of the rem- 
nant. We present multi-dimensional hydrodynamic simulations of the evolution of WD 
merger remnants under the action of an a-viscosity. We initializ e our c alculations us- 
ing the output of eight WD merger simulations from Dan et al. (2011), which span a 
range of mass ratios and total masses. We generically find that the merger remnants 
evolve towards spherical states on time-scales of hours, even though a significant frac- 
tion of the mass is initially rotationally supported. The viscous evolution unbinds 
only a very small amount of mass (< lCP 5 A/0). Viscous heating causes some of the 
systems we study with He WD secondaries to reach conditions of nearly dynamical 
burning. It is thus possible that the post-merger viscous phase triggers detonation of 
the He envelope in some WD mergers, potentially producing a Type la supernova via 
a double detonation scenario. Our calculations provide the proper initial conditions for 
studying the long-term thermal evolution of WD merger remnants. This is important 
for understanding WD mergers as progenitors of Type la supernovae, neutron stars, 
R Coronae Borealis stars and other phenomena. 

Key words: white dwarfs - hydrodynamics - supernovae: general 



> 

X 



1 INTRODUCTION 

Systems consisting of two white dwarfs (WDs) are natural 
outcomes of binary stellar evolution. These binaries are not 
static; absent any other torques the loss of angular momen- 
tum via gravitational wave (GW) emission will drive the bi- 
nary together. Programs such as the SWARMS survey publ 



can have a significant impact on the dynamics of the binary 



lally et al.||2009[ ) and the ELM survey ( |Brown et al.||2010[ ) 

have dramatically increased the number of known WD bi- 
naries, including some systems which will merge within a 

The Galactic population 



2012) 



Hubble time (Kilic et al 
of WD binaries is expected to be a source of unresolved 
GW foregrounds at mHz frequencies, though only a hand- 
ful of presently known systems would be individually de- 



tectable by a space-based GW interferometer mission ( Nclc- 
mans |[2u09| . 

Details of the inspiral, in particular whether tidal 
torques cause the binary to be synchronized and the loca- 
tion of the tidal heating, are active areas of inquiry that 



Einstein Fellow 



and the thermal state of the WDs (Fuller & Lai 20121. As 



the orbital separation shrinks, the less massive (and hence 
larger) WD will eventually overflow its Roche lobe and be- 
gin transferring mass to the companion. The stability of this 
mass transfer depends on e.g., whether the material forms 
a disc or flows directly onto the companion, which in turn 
depends on the mass ratio (q) and total mass (M to t) of the 
binary (e.g. |Marsh et al.|2004| ). 

Those systems which do undergo unstable mass transfer 
and subsequently merge have been of substantial theoreti- 
cal interest. In particular, such systems have received atten- 



tion as the possible progenitors of Type la supernovae ( Iben 
fc Tutukov||1984||Webbink||1984[ ). Considerable work exists 
exploring this "double degenerate" scenario and recent ob- 



servational results have begun to favor it (e.g. Bloom et al. 
2012 Schaefer & Pagnotta 2012 1. Another possibility is that 
double white dwarf binaries with total masses exceeding the 
Chandrasekhar mass undergo accretion induced collapse to 
form a neutron star (e.g. |Saio fc Nomoto 1985 1. Less massive 
double degenerate systems are likely to have non-explosive 
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outcomes and have been invoked to explain objects like the 



2 NUMERICAL METHODS 



R Coronae Borealis stars and extreme helium stars (Web- 
bink|p84l |Saio fc Jeffery|2000| |Clayton et aT1[2007l ). 



An accurate simulation of the merger process requires 
a 3D code without prescribed geometry and with good nu- 
merical conservation properties. For these reasons, the pi- 



oneering study of Benz et al. ( 1990 1 used smoothed par- 



ticle hydrodynamics (SPH). More recent studies (e.g. Dan 



et al. 2011 Raskin et al. 2012 Pakmor et al. 2012) have 



improved on these first results by contributing additional 
physics, more accurate initial conditions, higher resolution 
and more sophisticated numerical techniques. These simula- 
tions follow the evolution of the binary through the tidal dis- 
ruption of one of the components. In some cases the merger 



is sufficiently violent that an explosion may result ( Pakmor 
et al. 2010 Dan et al.||2012l). When the merger itself does 



not trigger an explosion, some material from the disrupted 
lower mass WD forms a shock-heated layer at the surface of 
the primary WD while the rest of the material forms a thick 
disc at larger radii. 

The evolution of such systems has frequently been 
treated in the literature as a long-lived (~ 10 5 yr) phase of 



accretion from a disc at the Eddington limit (e.g. Nomoto & 
Iben|1985| |. This picture was improved by |Yoon et al.| ( |2007[ ), 
who considered accretion at a similar rate but onto a hot en- 



velope, and by van Kerkwijk et al. (20101, who made simple 



a-disc estimates of the accretion time-scale and found it to 
be far more rapid (~ hours) than the time-scale for accretion 
at the Eddington limit. 



Recently, Shen et al. ( 2012 ) provided a new model of the 



different evolutionary phases of WD merger remnants. They 
argued that the evolution is much more "star-like" than the 
accretion disc oriented models that have dominated the lit- 



erature. More concretely, Shen et al. (20121 showed that the 



rapid dynamical evolution of the merger (~ 10 s) gives way 
to a longer lived viscous phase driven by magnetohydrody- 



namic instabilities (~ 10 



s) before the onset of a long 



(~ 10 4 yr) thermal phase. In contrast with previous work, 
this implies that the long term evolution of a white dwarf 
merger remnant is not determined by accretion, but rather 
by the internal redistribution of heat/momentum and the 
external cooling rate of the viscously heated, nearly shear- 
free remnant. 



In Shen et al. (20121, the viscous evolution was calcu- 



lated in ID using a 7-law equation of state. The goal of this 
work is to refine the understanding of the outcome of the 
viscous evolution of WD merger remnants using higher di- 
mensional numerical simulations. In addition, we consider a 



wider variety of WD+WD systems than Shen et al. (20121, 



who focused on roughly Chandrasekhar mass CO+CO merg- 
ers. 

In §2 we outline the numerical methods we use, includ- 
ing how we construct our initial conditions from simulations 



by Dan et al. (2011 1. In §3 we present the results of each of 



our calculations. §4 provides a discussion of the end states of 
the calculations. In §5 we state our conclusions and propose 
avenues for future work. In an Appendix, we show various 
test calculations that confirm the results we focus on in the 
main text. 



We perform our calculations using the ZEUS-MP2 (Hayes 



et al. |2006) code, a massively parallel implementation of the 



algorithms used in the ZEUS family of codes. These codes 
solve the fluid equations using finite differences on a stag- 
gered mesh. The internals of ZEUS are well-documented in 
the literature (for example, Stone & Norman 1992 1. While 
there have been other, more modern developments in astro- 
physical fluid codes, we chose to use ZEUS-MP2 because of 
its supported features (e.g. spherical coordinates, non-ideal 
equations of state) and because its structure allows for the 
relatively easy addition of new features. 

Our calculations are done in spherical coordinates, an- 
ticipating the evolution of the remnant to a quasi-spherical 
end state. In order to minimize the computational demands, 
we primarily perform 2.5D simulations, in which vector 
quantities can have a <f> component, but its value does not 
vary along the <f> direction. In general, we also assume reflec- 
tion symmetry about 9 = tt/2. In the Appendix, we briefly 
report additional calculations which confirmed the validity 
of these simplifications. 

Our typical computational domain is characterized by 
the grid spacing in the r and 8 directions and by the radius 
of the inner boundary. Unless otherwise specified, we adopt a 
logarithmic radial grid with N r = 64 points per decade. The 
angular grid is uniform from [0,7r/2] with No = 48 angular 
zones. These values give a grid in which individual cells are 
roughly equal in radial and 6 extent. We choose an inner 
radius such that only 0.1 per cent of the mass lies interior to 
that radius and then place the outer boundary at 10 4 ri nncr . 
We perform higher resolution simulations to confirm that 



our simulations are converged (see [ Al I 



We make several modifications to ZEUS-MP2 (based off 
of v2.12) in order to perform our calculations; we describe 
these modifications in the rest of this section. 



2.1 Non-ideal Equation of State 

We modify the code by the addition of the non-ideal 
Helmholtz equation of state (EoS), which provides an 
electron-positron EoS valid over a large range of physical 
conditions combined with the equations of state for an ideal 



gas of ions and for blackbody radiation (Timmes & Swesty 



2000). With the addition of the EoS routines, one small al- 



gorithmic change is made: as suggested in | Stone &i Nor-| 



( 1992 I, a predictor-corrector method is used to improve 



energy conservation during the calculation of the compres- 
sional heating term. 



2.2 Shear Viscosity 

In order to approximate the effects of magnetic stresses, we 
add shear stress terms to the hydrodynamic equations. That 
is, we are solving the equations 

Dtp + /£&=<} (1) 
pDtVi = -drP - + djTij (2) 
pD t (e/p) = -PdiVi + T t3 T,j/(pv) (3) 

where D t = dt + Vidi is the convective derivative and we ob- 
serve the usual Einstein summation conventions. The pres- 
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sure is denoted by P, and the mass and internal energy den- 
sities are represented by p and e respectively. The velocity 
vector is Vi. The anomalous stress tensor Tij is defined as 



Tij — pv (diVj + djVi) 



(4) 



where v is the dynamic viscosity. 

A very similar modification of the ZEUS-2D code was 
made by |Stone et al.| ( |1999[ ). We benefited from inspecting 
the source code that was used to perform the calculations 
reported in that work. We also used the results reported in 
| Stone et al.| ( |1999[ > as a reference against which to test our 
own implementation. 

The viscous terms are evaluated using an operator split 



method and are updated during the source step (Stone & 



Norman||1992 1. To ensure numerical stability, the time step 
At must be chosen to be less than At v i sc ~ min((Ar) 2 /i/), 
where the minimum is evaluated over the computational do- 
main. 

As a dimensionally motivated form for the dynamic vis- 
cosity coefficient, we adopt 



n fc (r) 



(5) 



where c s is the local sound speed and Qk is the Keple- 
rian angular velocity calculated using the mass enclosed at 
a given spherical radius. Portions of the merger remnant 
(see Fig. [T]) are unstable to the magneto-rotational insta- 
bility (M RI; |Balbus fc Hawleyl|1991| ) and the Tayler-Spruit 
dynamo ( |Tayler 1973 Spruit ]2002 1 . These processes may 
generate viscous stresses corresponding to a in the range 
1(T 4 - IP" 1 ; for order of magnitude estimates, see Shenl 
(20121. We adopt a fiducial value of a = 3 x 10" 



et al. 



though we confirm that the results of our calculations are 
not sensitive to this choice (see j ]A2[ ). 

As one moves to small r, (the origin being at the centre 
of the surviving WD; see {2A), c s and Ofc approach con- 
stant values. We are using a logarithmic grid, so Ar oc r 
and therefore Ai v isc oc r 2 . The timestep constraint imposed 
by the Courant-Friedrichs-Lewy (CFL) condition depends 
linearly on Ar, so AtcFL oc r. At sufficiently small radii, 
the viscous timestep becomes much less than the timestep 
required by the CFL condition. In practice, this occurs at 
a radius which is in our computational domain. In order 
to evolve the remnant over many viscous times, we apply 
the following ad hoc prescription. Within some radius r v we 
suppress the viscosity by a factor of 1 /r such that the ratio 
of tviac/tcFL remains constant. In order to make the cutoff 
smooth, the exact prescription is 



1/rf 



(6) 



where /3 = 4 and r v is approximately the half-mass radius 
of the inner WD. As shown in the Appendix, we have ver- 
ified that our results are insensitive to the details of this 
prescription. Physically, we would not expect the viscosity 
prescription in Equation [5] to hold as r — > because this 
region within the WD is in approximate solid body rotation 
and is not MHD unstable. 

Local numerical simulations of the MRI in accretion 
discs indicate that the azimuthal components of the stress- 
, are roughly a factor of 10 larger than the 



tensor, T rl p and 



Our default assumption then is to evolve with only these 
components in the stress tensor being non-zero. In §A3| we 
test the effect of including all components and find some 
quantitative, though not qualitative, differences between the 
two choices. 



2.3 Nuclear Burning 

In several of the model systems the temperatures reached 
are sufficiently high that the energy release from fusion is 
not negligible on the viscous time-scale. However, the con- 
ditions are not such that the local burning time-scale ever 
falls below AtcFL- Therefore, in order to minimize the com- 
putational load associated with calculating the burning, we 
implement an extremely simple 5 isotope nuclear network 
which is explicitly integrated at the hydrodynamic timestep. 
This captures the bulk of the energy release. 

The 5 species we track are He, C, O, Ne, Mg (these are 
the 5 isotopes present in the initial compositions). These 
isotopes are connected by 4 processes, the triple-a process 
and a capture on each of C, O, and Ne. The rates of these 
processes were taken from the JINA REACLIB database 



(Cyburt et al.|2010l. We neglect additional physics such as 



screening corrections because the burning primarily occurs 
at densities where such effects are unimportant. We refer to 
our own burning implementation as HeCONe. 

As a test of both our own implementation and the as- 
sumptions that motivate it, we also coupled the 13 isotope 
a-chain network aproxl3 to the code ( jTimmes et al.|[2000[ ). 
The results were virtually identical, confirming the validity 
of our approach. See |A6| for quantitative comparisons of the 
results of these tests. 

ZEUS-MP2 provides the ability to advect scalar quan- 
tities; we use this feature to track the mass fractions of the 
isotopes in our network. The algorithms do not guarantee 
the sum of the mass fractions remains equal to one after 
the advection step. Methods to restore this constraint in 
fluid codes have been reported in the literature (for exam- 
ple, [Plew^&^ullcr 1999 1. However, for the sake of simplic- 
ity, immediately before evaluating the energy release from 
nuclear burning, we enforce the constraint ~^Xi = 1 by 
appropriately adjusting the mass fraction of the most abun- 
dant isotope. 



2.4 Construction of Initial Conditions 

Our starting point is the results of SPH simulations of double 



white dwarfs performed by Dan et al. (20111. The simula- 



other components (e.g. Hawley et al. 1995 Stone et al. 1996 1 



tions were notable for their use of a more accurate initializa- 
tion of the binary system at the onset of mass transfer than 
had been used in previous work. We calculate the viscous 
evolution of the merger remnants formed in each of their 
eight production runs, the parameters of which are summa- 
rized in Table[l] Throughout the rest of this paper, when we 
refer to "initial conditions," this refers to the matter config- 
urations at the end of these SPH simulations. 

In order to map between the output of the SPH calcu- 
lations and the staggered mesh of ZEUS-MP2 we adopt the 
following procedure. In SPH, the value of a quantity / at a 
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ID 


M 2 


Mi 


A/tot 


<; 


^cnd 


C 2 


Ci 


PI 


0.2 


0.8 


1.0 


0.25 


95.0 


He 


CO 


P2 


0.3 


1.1 


1.1 


0.27 


62.0 


He 


CO 


P3 


0.5 


1.2 


1.7 


0.42 


35.0 


He 


ONeMg 


P4 


0.3 


0.6 


0.9 


0.50 


50.0 


He 


CO 


P5 


0.6 


0.9 


1.5 


0.67 


35.7 


CO 


CO 


P6 


0.2 


0.3 


0.5 


0.67 


30.0 


He 


He 


P7 


0.3 


0.4 


0.7 


0.75 


18.0 


He 


He 


PS 


0.9 


1.2 


2.1 


0.75 


30.0 


CO 


CO 



Table 1. A summary of the systems simulated by |Dan et al.| 
(2011). ID is their run number. M 2 is the mass (in Mq) of the 
secondary, the less massive of the two WD; Mi is the mass of 
the more massive primary. Mtot is the total mass of the system 
and q = Mi/M\ is the mass ratio. t en d is the duration of the 
SPH simulation in terms of the initial orbital period. Ci and C 2 
are the compositions of the primary and the secondary WDs. See 
Table 1 in |Dan et aT] ( |2011^ for more details. 



position x is given by 



(7) 



where W is the kernel function. The quantity of interest 
associated with the i-th SPH particle is denoted The SPH 
particle has mass, density, position and smoothing length 
rrii, pi, Xi and hi respectively (e.g. Monaghan"| [l992 ) . The 



sum runs over the total number of SPH particles, n. 

Schematically, we construct our grid-based initial con- 
ditions by evaluating the five quantities that ZEUS-MP2 
evolves (mass density p, internal energy density e and veloc- 
ity v) at each grid point. In our standard 2D simulations, 
we construct initial conditions that are explicitly invariant in 
the 4> direction and are reflection-symmetric about 9 — tt/2. 
Given these conditions, the total linear momentum of the 
remnant is guaranteed to be zero. Therefore, we choose the 
origin of our simulation coordinate system to be at the point 
of peak density, corresponding to the centre of the more mas- 
sive (surviving) WD. 

Explicitly, in order to calculate the value of a density (e 
or p) at a grid point with coordinates (ri,9j) we evaluate 



Pij 



2JV, 



jrr [PSPH(ri, 6j, <j>k) + pspn(n, tt ~9j,<f) k )] (8) 



where N$ is the number (typically N$ = 32) of equally- 
spaced points used to cover the interval <j> £ [0, 2-7r). Now 
and throughout this section, the subscript SPH indicates a 
quantity extracted from the SPH simulation by the evalua- 
tion of Eq. [7] 

Constructing the initial velocity vector requires slightly 
more complicated <j> averaging. From the SPH data, we first 
construct the full velocity vector (with the components rep- 
resented in Cartesian coordinates) at each of the grid loca- 
tions where a component of the velocity will be defined. The 
staggered mesh employed by ZEUS-MP2 means that each 
velocity component is defined at a different spatial point 



(for details see Stone & Norman 19921. We take this into 



account, though we do not manifestly indicate it in the for- 
mulae below for the sake of compactness. 

It is not simultaneously possible to conserve the kinetic 
energy of the fluid and its linear and angular momentum 



when performing the (^-averaging. (This is simply a state- 
ment of the fact that in a non-uniform field, (u 2 ) 7^ (v) 2 .) 
Given that we are interested in investigating the angular 
momentum evolution of the remnant, we choose to conserve 
momentum. In practice, the difference is relatively small; for 
the fiducial remnant, this averaging procedure changes the 
total kinetic energy by 1 per cent. 

Therefore, to obtain the value of a component of the 
velocity v, defined at a point (r»,0j) we calculate 

1 1 N * 

ZJV Pn fe= i l y J 



+ PSPH(ri,7T - 



,4>k) ■ e»j(0k)] 



where p is the momentum vector and e%j is the unit vector 
of the velocity component of interest at the point. 



3 RESULTS 

We expect that systems with similar mass ratios (q) and 
total masses (Mtot) will undergo similar evolution. Since 
the composition of a WD maps to a relatively well-defined 
mass range, we organize our discussion by the composition 
of the merging WDs. First, we discuss our fiducial 0.6+0.9 
M Q CO+CO system. The outcomes of He+He, He+CO and 
He+ONeMg mergers are discussed on a more limited basis, 
emphasizing only the notable differences between these sys- 
tems and our fiducial one. See Table [2] for a summary of the 
properties of our primary simulations. As shown in Table [l] 
Dan et al. (2011 1 label their simulations with a short iden- 



tifier of the form Pn, where P represents production and n 
in an integer. For our own short IDs (shown in Table ^J), 
we simply prepend Z (representing "ZEUS") to the ID of 
the Dan et al. (2011) simulation that was used as the initial 



conditions. 

In order to easily visualize our multi-D simulations, we 
will plot spherically averaged quantities against spherical 
radius. To calculate densities we perform a volume average, 
e.g. 

1 



p(r) = 



dO sm(9)p(r,9) 



(10) 



so that the appropriate quantity (e.g. mass) is conserved. 
To calculate other thermodynamic quantities such as tem- 
perature, we first calculate the spherically averaged mass 
and energy densities and then apply the equation of state. 
To calculate angular velocities, we restrict the average to a 
45° wedge centred on the equator. As our simulations evolve 
toward a spherical end state, these ID averages become an 
increasingly complete summary of the 2D structure of the 
remnant. 



3.1 CO+CO systems 

Our fiducial system (ZP5) is a super-Chandrasekhar 
CO+CO merger with Mtot = 1.5M and q = 2/3. We simu- 
late this system for 3 x 10 4 s, which is ~ 5 x 10 7 timesteps of 
the hydrodynamics code. The simulation conserves mass to 
1 part in 10 4 , energy to 0.5 per cent and angular momentum 
to one part in 10 3 . The evolution of an identical system was 



discussed in Shen et al. (2012), who performed a simple ID 
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ID 


M 2 -4 


- Ml 


^"innc 


[cm] 


Network 


tend [s] 


ZP1 


0.2 + 0.8 


3.6 


xlO 7 


HeCONe 


3.0 xlO 4 


ZP2 


0.3 - 


b 1-1 


2.2 


xlO 7 


HeCONe 


3.0 xlO 4 


ZP3 


0.5 - 


h 1.2 


1.9 


xlO 7 


HeCONe 


1.5 xlO 4 


ZP4 


0.3 - 


h 0.6 


4.4 


xlO 7 


HeCONe 


4.0 xlO 4 


ZP5* 


0.6 - 


b 0.9 


2.9 


xlO 7 


None 


3.0 xlO 4 


ZP6 


0.2 - 


1- 0.3 


6.6 


xlO 7 


HeCONe 


4.0 xlO 4 


ZP7 


0.3 + 0.4 


5.7 xlO 7 


HeCONe 


4.0 xlO 4 


ZP8+ 


0.9 - 


h 1.2 


1.9 


xlO 7 


aproxl3 


1.0 xlO 4 



Table 2. Details of the viscous evolution calculations discussed in 
this paper. As an ID, we simply prepend Z (representing "ZEUS") 
to the ID of the |Dan et al.| ( |2011| > simulation that was used as 
the initial conditions. M2 + Mi is the mass of the secondary + 
primary in Mq . We will sometimes refer to systems by this sum. 
dinner is the location of the inner boundary of our computational 
domain. Network indicates which nuclear network was used in 
the calculation. t en( j is the end time of the simulation. All of the 
simulation parameters which were held fixed across all runs are 
discussed in the text. * fiducial model discussed in the most detail 




in the main text (in f 3.1 I t this simulation had a lower resolution, 
N r ,N g = 48,32. 



M r [M G 



calculation of the viscous evolution. Our multidimensional 
calculations confirm the schematic picture presented therein. 

At the end of the SPH simulation, the primary white 
dwarf is relatively undisturbed and is surrounded by the 
remnants of the disrupted secondary. More than half of the 
disrupted material has primarily rotational support; the re- 
mainder was shock-heated in the merger and has thermal 
support. (A small amount ~ 10 Mq is unbound in a tidal 
tail.) The merger remnant is in quasi-hydrostatic equilib- 
rium, which we confirm by evolving these initial conditions 
without the action of viscous torques for many dynamical 
times, observing little change. 

The black lines in Fig.[l]show the initial rotation profile. 
The primary WD is rotating rigidly; exterior to that is the 
disrupted secondary which is in Keplerian rotation. This hot, 
fully ionized material is unstable to magnetohydrodynamic 
instabilities such as the MRI. The turbulence generated by 
the saturation of the MRI leads to an enhanced viscosity 
and concomitant transport of angular momentum to larger 
radii ( [Balbus fc Hawley|199lT|Balbus|2003| ). As described in 
§2.2| we model this using an a- viscosity. 

The viscosity also liberates energy present in the <j>- 
velocity shear. Fig.[2]shows the evolution of the temperature 
and specific entropy profiles during the viscous evolution. 
Fig. [3] shows the final p — T distribution and the evolution 
of the temperature peak in the p— T plane during the viscous 
evolution. The viscous heating causes the peak temperature 
in the remnant to increase over the duration of the viscous 
phase by a factor of ~ 2, to ~ 8 x 10 8 K. The tempera- 
ture peak is at a density of ~ 5 x 10 5 gcm -3 and at those 
conditions the energy released from carbon burning exceeds 
neutrino losses and the burning becomes self-sustaining. 

The carbon burning in our fiducial model is an unimpor- 
tant energy source on the viscous time-scale, so the viscous 
evolution is not directly affected. However, the fact that car- 
bon burning begins during the viscous evolution means that 
a convective carbon burning shell will develop in ~ 10 6 s. 
One consequence of this is that we expect the material exte- 
rior to the temperature peak at the end of the viscous evolu- 



Figure 1. The evolution of the rotation profile of the fiducial 
0.6+0.9 Mq CO+CO remnant. The solid lines are the angular 
velocity Q and the dashed lines are its ratio to the Keplerian 
angular velocity, Q/Q^. (The angular velocity is calculated using 
the material in a 45° wedge centred on the equator.) In the initial 
profiles (black), most of the material from the disrupted secondary 
(M r > 0.9) is rotationally-supported with an angular velocity 
profile unstable to the MRI. The action of viscosity drives more of 
the remnant to solid body rotation and the accompanying heating 
leads to more of the remnant being thermally supported. We set 
a = 0.03 for all the results in the main text. The Appendix shows 
that our results are nearly independent of a. 



tion to quickly form a convective envelope. Future work will 
investigate the structure of this envelope, which is important 
for understanding the thermal evolution of the remnant and 
characterizing its observational properties. 

One of the most striking results of our multi-D simu- 
lation is that the merger remnant evolves towards a final 
quasi-spherical steady state. (Given that there is rotation, 
the final state will actually be oblate, though in practice, 
the rotational velocities we find imply that it is quite spher- 
ical.) To quantify this, we define a simple "aspect ratio" as 
follows: draw an isodensity contour and take the ratio of the 
distance from the origin at the equator to that at the pole. 
As a rule of thumb, we find that the aspect ratio associated 
with a given radius approaches unity after about 10 viscous 
times have passed at that radius. The bottom panel of Fig. 
[4] shows this convergence clearly. The primary WD ends up 
with a thermally supported, nearly spherical envelope. The 
top panel of Fig. [4] shows the mass enclosed as a function of 
radius which illustrates how the outer thermally supported 
envelope expands to larger radii during the viscous evolu- 
tion. 

Our multi-dimensional simulations also allow us to cap- 
ture processes like convection. We find that the viscous heat- 
ing causes the remnant to evolve towards a convectively un- 



stable state. Recently, Garcfa-Berro et al. (20121 discussed 
how convectively generated magnetic fields in merger rem- 
nants could potentially explain the origin of high-field WDs. 
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Figure 2. The evolution of the temperature (top) and specific en- 
tropy (bottom) profiles of the fiducial 0.6+0.9 M CO+CO rem- 
nant. Viscous heating increases the peak temperature by roughly 
a factor of two. Convection and viscous heating both contribute 
to the entropy evolution of the material from the disrupted sec- 
ondary. As discussed in the text, these curves are ID spherical 
averages of our 2D simulations. 
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Figure 3. The evolution of the temperature peak in the p — T 
plane. The dotted line indicates the break-even point where the 
energy release from carbon burning is equal to neutrino losses. 
The filled square (circle) is the peak temperature and correspond- 
ing density at the start (end) of the simulation, and the dashed 
line that connects them traces its evolution. The solid line is the 
full ID p — T profile of the quasi-spherical end state. The grey 
dash-dot line indicates where gas and radiation pressure are equal. 
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Figure 4. The evolution of the mass enclosed (top) as a func- 
tion of spherical radius and the aspect ratio (bottom; see text for 
definition) for the fiducial 0.6+0.9 Mq CO+CO remnant. The 
convergence of the aspect ratio to a constant value PS 1 provides 
an indication that the end state has approximately been reached. 

They noted that the conditions at the end of their own 
SPH simulations were unstable by the Schwarzchild crite- 
rion. However, given that this system has substantial rota- 
tional support and that we are evolving it in axisymmetry, 
a more appropriate test is the H0iland criterion. Our initial 
conditions are not unstable by the H0iland criterion and do 
not to evolve towards an unstable state without the action 
of the viscous stresses [£] 

Fig. [5] shows the 2D evolution of the fiducial system. In 
addition to the entropy and temperature, we plot two energy 
densities which are helpful in interpreting the evolution. One 
is the free energy available in the 0-velocity shear 

KE^^ip^-^) 2 , (11) 

which shows the energy available for viscous heating. The 
other is the kinetic energy density in non-azimuthal motions 

KE r , fl = \p [y\ + v 2 e ) , (12) 

which is related to the kinetic energy associated with con- 
vective motions. Fig. [5] also plots isodensity contours, which 
emphasize the approach of the remnant to a spherical state. 

Given previous work on viscous, geometrically thick ac- 
cretion flows, one might expect that material would out- 
flow along the poles during the viscous evolution. When the 
viscous time is much less than the cooling time and the 
mass inflow is assumed to be conservative (that is, mass 
does not leave the system), the transport of energy and 
the release of gravitational potential energy are such that 

1 The physical initial conditions are of course unstable in MHD, 
as it is the MRI that is giving rise to the effective viscosity. 
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Figure 5. A visual overview of the 2D evolution of the fiducial 0.6+0.9 Mq CO+CO merger remnant. Within each panel, the top 
two subpanels are thermodynamic quantities (s, T) and the bottom two subpanels are kinetic energy densities (non-azimuthal, </>-shear). 
The black contours are density, spaced one per decade. The dashed contour is p = 10 3 g cm -3 . The main panels are snapshots of the 
simulation at the indicated times. Top Panel: The initial conditions, note the large "free" energy apparent in the shearing, Keplerian disc. 
Middle Panel: The action of viscosity has dissipated some of the shear and heated the material. The remnant has become convectivcly 
unstable as can be seen in the striation of the non-azimuthal KE. Bottom Panel: The remnant has settled down into a quasi-spherical 
steady state. 
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material has a positive Bernoulli parameter (Be) (Narayan 



& Yi 1994 Blandford & Begelman 19991. Therefore, solu- 



tions in which the mass flow is not conservative may be 
more physical. Non-radiative accretion flows are also pre- 
dicted to be convectively unstable. Models based on these 
ideas (e.g |Blandford fc Begelman|2004| | developed solutions 
with prominent outflows. Hydrodynamic numerical simula- 



tions such as those by Stone et al. ( 1999 1 exhibited the slow 



outflow of marginally bound material in the polar direction. 



MHD simulations such as those by Stone & Pringle (20011 



found somewhat more prominent outflows than in the hy- 
drodynamic simulations. 

We do not observe outflows in our simulations. Fig. 
[6] shows the fraction of mass on our grid with positive 
Be = /Be>o- Initially /boo ~ 3 x 10~ 3 , corresponding to 
the unbound material in the tidal tail. This material flows 
out of the outflow boundary and afterwards there is little 
unbound mass (/boo ~ 10 ). In order to isolate the ef- 
fects of viscosity on the unbound material, we ran a sim- 
ulation without the viscosity and calculated /boo (dotted 
blue line). The orange line in Fig. [fj] shows the integrated 
difference in the mass that flowed through the outer bound- 
ary with Be > in simulations with and without viscosity. 
This difference is very small < 10 _5 Mq. We do not claim 
that this specific value is robust, but the conclusion that the 
viscous evolution of WD merger remnants unbinds a only a 
very small amount of mass appears to be. 

The initial conditions of our simulations are rather dif- 
ferent than the initial conditions of most radiatively ineffi- 
cient accretion simulations. Such simulations typically allow 
the material to move through several orders of magnitude 
in radius before reaching an inflow boundary representing a 
black hole. By contrast, the radial dynamic range between 
the surface of the primary WD and the bulk of the material 
in the initial rotationally-supported disc is small, a factor 
of ~ 5. The presence of a "hard surface" (the primary WD) 
means that as material accretes, the radius where material is 
pressure-supported moves outward, further suppressing the 
dynamic range between the effective inner boundary and the 
disc. 

In order to understand the results of our WD merger 
remnant simulations in the context of accretion tori simu- 
lations, we generate accretion tori like those in |Stone et al.| 
( 1999 1 and adjusted the inner boundary condition (reflect- 
ing vs. inflow) and the dynamic range between the initial 
torus and the inner boundary. We run the simulations for 
several orbits and calculate the resulting amount of unbound 
material /boo- At a radial dynamic range between the ini- 
tial torus and the inner boundary of 100 and with an inflow 
boundary, we find /boo ~ 4 x 10~ 2 , much larger than in 
our WD merger remnant calculations (see Fig. [6j. Decreas- 
ing the dynamic range to 10 results in /boo ~4x 10 -3 . At 
this dynamic range, changing the inner boundary condition 
to reflecting causes /boo 1° peak ~ 10 -4 and then fall as the 
simulation continues. These results support our conclusion 
that only a very small amount of mass is unbound during the 
viscous evolution of WD merger remnants (/boo 10 ) . 

In addition to our fiducial 0.6+0.9 M Q CO+CO sys- 
tem, we also simulate a very super-Chandra 0.9+1.2 Mq 
system (ZP8). This system quickly starts C+C burning, al- 
though the burning does not become dynamical (see Section 
4.2 1. While the energy release from the burning on the vis- 
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Figure 6. Viscously unbound material. The left (blue) scale in- 
dicates the fraction of mass with positive Bernoulli constant at a 
given time in the evolution of the fiducial 0.6+0.9 Mq CO+CO 
merger remnant. The dominant contribution is the tidal tail; the 
large decrease in /boo over the first 5000 s is this material flow- 
ing out of the simulation domain. The solid (dotted) line is the 
mass fraction with Be > in a simulation with (without) vis- 
cosity. The right (orange) scale is the integrated amount of mass 
that has flowed out of the simulation domain with positive Be due 
to the influence of viscosity. This is the integrated difference be- 
tween the two blue curves. Little additional material (^ 10 Mq) 
is unbound by the viscous evolution . 



cous time is locally non-negligible, the mass outflow is not 
affected by the presence of nuclear energy generation. As 
shown in the Appendix, the energy release is not significant 
enough to affect the global behavior of the remnant. 



3.2 He+He systems 

The evolution of He+He merger remnants is broadly similar 
to our fiducial CO+CO case. The larger size and lower mass 
of the He WDs mean that the temperatures reached at the 
end of the SPH calculations are not as high. However, these 
temperatures are still high enough that we elect to track the 
energy release from He burning in our simulations. We do 
this using the simple nuclear network described in Section 
1231 

Fig. [7] shows the evolution of the temperature and ro- 
tation rate for a 0.2+0.3 Mq system (ZP6). This has the 
same mass ratio q = 2/3 as the fiducial system, but with a 
total mass 3 times lower. The initial temperature and rota- 
tion profiles look qualitatively similar to our fiducial system. 
Appropriately scaling these values by the total mass, they 
are even quantitatively similar. However, the final state ap- 
pears somewhat different, most conspicuously because of the 
narrow temperature peak that forms at an enclosed mass of 
~ O.38M . 

This qualitative difference in evolution is most clear in 
Fig. [8] which shows the evolution of the temperature maxi- 
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Figure 7. A ID summary of the 0.2+0.3 Mq He+He remnant 
evolution. The temperature (top panel) and rotation (bottom 
panel) profile at beginning, intermediate and final times. The 
temperature evolution appears qualitatively different than our 
fiducial model; as explained in the text, this is due to this lower 
mass remnant remaining gas pressure dominated. The rotation 
evolution is qualitatively the same as the fiducial model. 

mum and the corresponding density. The temperature max- 
imum evolves to lower density during the viscous phase, un- 
like in CO+CO mergers where it evolves to higher density 
(see Fig. j3j). This effect is unrelated to the presence of fu- 
sion, as the time-scale for burning is still relatively long. 
This difference in evolution is also not a qualitative differ- 
ence in the merger outcome or the viscous evolution, but 
rather is due to the different contribution of gas and radi- 
ation pressure in the merger remnants. The grey dash dot 
line in Figs. [3] and [8] shows where P gas = P ra d- The He+He 
case remains gas pressure dominated, so the location of the 
peak temperature corresponds to the location of the vis- 
cous heating. Lower densities, which are at larger radii and 
have corresponding longer viscous time-scales, are heated at 
later times. In the CO+CO case, radiation pressure has a 
larger relative contribution initially and once material has 
been heated such that radiation pressure begins to become 
important, additional heating is no longer as effective in rais- 
ing the temperature. Therefore, larger relative increases in 
T occur at higher densities where gas pressure continues to 
dominate, and thus the peak T moves to higher densities. 




1 10 10 2 10 3 10 4 10 5 10 6 
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Figure 8. The evolution of the temperature peak for the Hc+Hc 
merger remnants. The dotted line labeled t^ UIn jj c = 10 4 s indi- 
cates the region where the time-scale for energy release by He fu- 
sion is equal to the time-scale of the simulation. The filled square 
(circle) is the peak temperature and corresponding density at the 
start (end) of the simulation, and the dashed line that connects 
them traces its evolution. The solid line is the full ID p — T profile 
of the quasi-spherical end state. The grey dash-dot line indicates 
where gas and radiation pressure are equal. 



The 0.3+0.6 M He+CO system (ZP4) is the only one 
of the systems we simulate in which the final state devi- 
ates significantly from approximate spherical symmetry. The 
remnant itself is spherical, but at the interface between the 
material from the He and CO WDs the composition varies 
between the equator and pole. This explains why the fi- 
nal peak temperature does not lie on the final spherically- 
averaged p — T profile shown in Fig. [9] This system has the 
highest mass ratio of any mixed composition system we sim- 
ulated and at higher mass ratios the secondary more strongly 
affects the primary. However, because the efficiency of mix- 
ing likely depends on dimensionality and angular momentum 
transport (a- viscosity vs MHD), we do not expect our work 
to provide a robust prediction of the details of such mixing. 



3.3 He+CO systems 



3.4 He+ONeMg 



In the He+CO systems, the primary is more massive and 
hence more compact than in the He+He mergers. This leads 
to higher temperatures during the merger. The lower tem- 
peratures required for He burning (versus C burning) mean 
that the effects of nuclear burning are more pronounced in 
these systems. Fig.[9]shows the evolution of the temperature 
peak in these simulations. 



One system in our study is composed of a 0.5Mq He WD 
and a 1.2M© ONeMg WD. Its evolution is very similar to 
the systems previously discussed. Notably, the high primary 
WD mass means that the He (from the secondary) reaches 
quite high temperatures. The burning time in this system is 
thus quite short, less than the viscous time, though not less 



than the dynamical time (see [4.21 
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Figure 9. The evolution of the temperature peak for the He+CO 
remnants. The dotted line labeled erjC = e„ indicates the break- 
even point where the energy release from carbon burning is equal 
to neutrino losses. The dotted line labeled t^ UIn h c = 10 4 s indi- 
cates the region where the time-scale for energy release by He fu- 
sion is equal to the time-scale of the simulation. The filled square 
(circle) is the peak temperature and corresponding density at the 
start (end) of the simulation, and the dashed line that connects 
them traces its evolution. The solid line is the full ID p — T 
profile of the quasi-spherical end state. The red dotted section 
shows where the helium mass fraction exceeds 50 per cent. The 
grey dash-dot line indicates where gas and radiation pressure are 
equal. The peak temperature in the 0.3+0.6 p — T profile does 
not correspond to the final peak temperature indicated by the 
solid circle. At the spherical radius of the temperature peak, the 
chemical composition varies from pole to equator and hence the 
averaged temperature at that point is not equal to the 2D peak. 

4 DISCUSSION 
4.1 Fitting Formulas 

The end states of our simulations will be useful as a starting 
point for future work concerning the thermal evolution of 
WD merger remnants. To aid such work, we provide fitting 
formulae that allow one to easily construct a physical state 
which is in rough quantitative agreement with our results. 

The ID profiles we extract at the end of our calculations 
have the following schematic form. At the centre is a core 
of cold, degenerate material. This is surrounded by a hot 
envelope, the outer portion of which was convective and so 
has an entropy that is roughly spatially constant. 

This picture allows a simple, post hoc model of the end 
state of our simulations. We write down a piecewise equation 
of state in which there is a central mass (M c ) described by 
a zero-temperature equation of state. This is surrounded by 
an isothermal region corresponding to the temperature peak 
which has mass M tp . The rest of the external material has 
a polytropic equation of state. Empirically the polytropic 
index of n = 3 provides a good fit to all of our simulations. 
For systems at high temperatures, such as our 0.3+1.1 Mq 



ID 


pc [gem 3 ] 


M c 


M tp 


ZPlt 


8.8xl0 6 


0.71 


0.10 


ZP2t 


4.7xl0 7 


0.98 


0.12 


ZP3t 


9.5xl0 7 


1.05 


0.16 


ZP4t 


3.8xl0 6 


0.53 


0.13 


ZP5 


2.8xl0 7 


0.84 


0.20 


ZP6 


6.4xl0 5 


0.28 


0.08 


ZP7 


1.5xl0 6 


0.38 


0.12 


ZP8 


3.3x10 s 


1.11 


0.24 



Table 3. The parameters from our fits (see Equations 1 13| 1 5] and 
surrounding discussion). ID is the run ID. p c is the central density 
extracted from the end of our simulations. M c is the amount of 
mass (in Mq) in the zero-temperature degenerate core. Mt p is 
the amount of mass (in Mq) in the isothermal region, loosely 
corresponding to the temperature peak. ^marks those systems 
which have a mixed chemical composition. 

He+CO merger, this is unsurprising as the material in the 
convective region is nearly radiation dominated, implying 
an adiabatic index near 7 = 4/3. For systems such as low 
total mass He+He mergers (ZP6 & ZP7), the matter is gas 
pressure dominated, which would imply an adiabatic index 
of 5/3. However these systems have larger residual entropy 
gradients, such that the relationship P oc p 4//3 roughly holds. 
Since we can provide satisfactory fits without introducing an 
additional parameter, we choose n = 3 for all our fits. 
Quantitatively 

[Pzt(p) if M r < M c 
P(p) = I Kip if M c < M r < M c + Mtp (13) 

[aV +1/ " ifM c + M t p<Mr 

where Pzt is the pressure of a zero temperature Fermi gas 
with p e = 2 (e.g |Shapiro fc Teukolsl^[T983| ). Ki and K 2 
are set by the condition that p, P are continuous at the 
transitions between regimes. 

In combination with the equations of hydrostatic equi- 
librium and mass conservation in ID spherical coordinates 

dM r ,9 / 1 a \ 

~ = 47rr p (14) 
dr 

dP _ GM rP 

-fr - -5- ( lb > 

and a central boundary condition p(ri nnor ) = p c , this is suf- 
ficient to fully specify a ID model. We set p c to be the value 
of the central density at the end of our simulations. 

For each of our simulations we fit for the two masses 
M c and Mtp. Table [3] reports the results of these fits. Fig. 
[lO] shows the results of the fit to our fiducial model. The fit 
reproduces the observed quantities to within ~ 30 per cent. 
The fit is worst in the region described by the isothermal 
equation of state, which is unsurprising since this is the least 
physically motivated part of our effective equation of state. 

Our fitting procedure does not use or provide any spa- 
tial information about the chemical composition. As a rough 
approximation, one can simply retain the initial Lagrangian 
composition of the system with the secondary outside of the 
primary. In the mergers where the chemical compositions 
of the WDs were initially identical, this is a good approxi- 
mation because nuclear reactions do not significantly alter 
the composition (for the set of mergers we considered). For 
mergers where the WDs had different compositions (which 
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Figure 10. The simple two-parameter fit to the fiducial 0.6+0.9 
Mq CO+CO merger remnant (see Equat ions | 1 3) 1 5| and surround- 
ing discussion). The upper two panels show the pressure and den- 
sity profiles from the simulation as solid lines. The fits are shown 
as dashed lines. The bottom panel shows the relative error be- 
tween the fit and the simulation. The vertical grey lines show the 
position of the transitions in the piecewise equation of state. 



are marked in Table pjl, the assumption that the composi- 
tion is conserved in a Lagrangian sense is substantially more 
crude because of mixing and the effects of nuclear burn- 
ing. In those cases, these simple fits would be inappropriate 
for work in which inaccuracies in the chemical composition 
could have a large effect. 



4.2 Burning Time 

Recently there has been considerable interest in the possibil- 
ity of central carbon detonations triggered by the detonation 
of a helium layer on the surface of a CO WD. During a WD 
merger, conditions for detonations may be reached in in- 
stabilities in the accretion stream (jGuillochon et aL||2010[ ) 



or at surface contact (e.g. Dan et al. 20121. The systems 



we consider did not reach detonation conditions during the 
SPH simulations (though those could not resolve accretion 
stream instabilities). 

During the phase of evolution that we simulate, viscous 
heating does increase the temperature and either initiate or 
increase the rate of burning. Fig. [Tl] shows the minimum 
burning times and corresponding temperatures for the eight 
systems we simulate. We calculate the minimum nuclear 
burning time as tbum = cpT/e nuc , where cp is the specific 
heat at constant pressure, T is the temperature and e nuc is 
the specific energy generation rate from nuclear reactions. 

The minimum burning time is not necessarily located 
at the location of peak temperature, as differences in the 
chemical composition (for example, the presence of helium) 
may make the rate of energy release greater at a different 
location. In general, the viscous heating causes a monotonic 



increase in the temperature. Therefore, initially the burning 
time drops. Then, in cases where the burning time is less 
than the viscous time, changes in the composition (such as 
the depletion of helium) begin to shift the minimum burn- 
ing time to slightly lower temperatures where more material 
remains to burn. 

In the case of the 0.5+1.2 Mq merger, the burning time, 
which is ~ 40tdyn at the beginning of the simulation, de- 
creases to as low as ~ 2td yn , where the dynamical time is 
calculated as t^yn = P/(pgc s ). The ratio tbum/tdyn S; 1 
provides a rough criterion for possible detonation. Detailed 
conditions for detonations are still a topic of current research 
and almost certainly require resolving smaller length scales 
than our current simulations can do (for example, see dis- 
cussion m §3.2 of |Woosley fc Kasen|201lj ). 

Using the value of tbum/^dyn as a guide, viscous heating 
does not cause any of the remnants that we simulated to ex- 
perience dynamical burning. However, with the low value of 
ibum/idyn f° r the 0.5+1.2 Mq system and the temperature 
sensitivity of nuclear reactions, we expect that systems only 
slightly more massive would experience dynamical burning. 
Furthermore, if there are stochastic fluctuations, it is even 
possible that this particular system could experience dynam- 
ical burning. Several figures in Dan et al. ( 2012[ ) draw con- 
tours of £bum/idyn in the space of primary /secondary WD 
mass. For systems where the merger has brought the sys- 
tem close to the critical value of tbum/idyn, we anticipate 
that viscous evolution would cause this value to fall below 
1. Viscous evolution reduces tbum/idyn by a factor of ~ 20 
in run ZP3, so we speculate that dynamical burning might 
occur for those systems with 1 < tburn/*dyn 10 in figure 



8 of Dan et al. (2012 1. Future work will quantitatively ad- 



dress this question by simulating a wider range of merger 
remnants. 



5 CONCLUSIONS 

The merger remnants of binary white dwarfs are differen- 
tially rotating and unstable to MHD instabilities like the 
MRI. As outlined by |Shen et aL] (|2012[), MHD stresses 



give rise to a viscous phase of evolution which occurs on a 
time-scale much less than the thermal time. To investigate 
the outcome of this viscous evolution, we perform multi- 
dimensional hydrodynamic calculations of the evolution of 
WD binary remnants under the action of an a- viscosity. The 
initial conditions for these calculations are the SPH simu- 



lations by Dan et al. (20111. We find that these remnants 



evolve towards spherical states on time-scales of hours. This 



confirms the arguments in Shen et al. (20121 that the post- 



merger evolution of WD merger remnants is via viscous re- 
distribution of angular momentum that leads to nearly solid 
body rotation. The transport of angular momentum out- 
wards removes rotational support from the majority of the 
mass leading to a nearly spherical remnant. The dynamics 
during this phase is not consistent with accretion at the Ed- 
dington limit, as in previous models of WD merger remnants 



(e.g.|Nomoto fc Iben|1985l|Saio fc Nomoto|1998||2004||Pier 



santi et al. 2003b|a l. Instead, the viscous evolution of WD 
merger remnants is much more analogous to that of a dif- 
ferentially rotating star. 

Viscous heating associated with the approach to solid 
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Figure 11. The shortest burning time (top panel) and corre- 
sponding temperature (bottom panel) for each of our simulated 
systems. The x-axis is the mass of the primary WD. Two systems 
have the same primary mass of 1.2 Mq and are slightly offset 
on the x-axis for visual clarity (ZP3 to the left and ZP8 to the 
right). In the top panel, the circle represents the shortest burning 
time reached overall, that is at any point during the simulation; 
the cross represents the burning time at the end of the simulation. 
They are connected by a dashed line to guide the eye and indicate 
that intermediate values are achieved. The same symbols in the 
bottom panel show the temperatures at the corresponding loca- 
tions. Because of varying chemical composition, the temperature 
associated with the shortest burning time is not necessarily the 
global peak temperature. The right axis of the top panel and the 
orange circles show the ratio tbumAdyn ( as defined in the text) 
at conditions corresponding to the black circles. In no case does 
the burning time ever reach the dynamical time, though in ZP3 
it is within a factor of two. 



body rotation unbinds only a very small amount of mass 
(< 1(T 5 M Q in our fiducial calculation) . This is in contrast to 
some of the intuition developed in the context of radiatively 
inefficient accretion flows, which predict outflows. To un- 
derstand this, we perform simple accretion tori calculations 
which indicate that the relatively small radius difference be- 
tween the disc and the surface of the WD can explain why 



only a small amount of mass becomes unbound (see 6 3.1 1. 

Viscous heating causes one of the systems we simulate 
to reach conditions of nearly dynamical He burning, so it 
is possible that the post-merger viscous evolution triggers a 



detonation in some cases. Recently Dan et al. (20121 pre- 



sented a suite of more than 200 WD merger simulations 
which more thoroughly populate the q-Mt ot plane. They 
found that many of these systems reached the conditions for 
detonation during the merger (see for example their figures 6 
& 8). In our calculations, min(iburn/idyn) decreases by a fac- 
tor of ~ 10 during the viscous phase (see §4.2). As a result, 
we speculate that systems that reach tbum/idyn < 10 dur- 
ing the merger may reach conditions for detonation during 
the subsequent viscous phase. Such a system is not among 



those that we simulate (though run ZP3 was close). It is 
thus likely that heating during the post-merger viscous evo- 
lution significantly increases the region of parameter space 
in which one expects explosive burning during WD mergers, 
which could increase the rates of Type la or .la supernovae 
( jBildsten et al.|2007 l associated with WD mergers. 

Our purely hydrodynamic simulations cannot address 
the effects of magnetic fields. MHD simulations resolving 
the action of the MRI would allow a more realistic treat- 
ment of the viscous stresses than an a-viscosit^ though 
the quantitative insensitivity of our results to the value of a 
leads us to think that our conclusions are robust. Convert- 
ing our fiducial value of a to a magnetic field strength gives 
\B\ ~ \J Airapcg ~ 10 10 G. The implications of this estimate 
for the subsequent evolution of the merger remnant depend 
on the structure of the field. The generation of a large-scale 
field could lead to the formation of a strongly magnetized 
WD, which would be rapidly rotating and would quickly 
spin down via a magnetized wind. The presence of a strong 
magnetic field would also affect the conduction of heat in 
the interior of the WD. Alternatively, it is possible that the 
strong field is relatively small scale and so efficiently redis- 
tributes angular momentum in the interior of the remnant 
but does not significantly affect its global properties. 

The end states of our calculations provide a starting 
point for investigations of the long-term thermal evolution 
of WD merger remnants. In our fiducial case, we expect that 
the luminosity from the nuclear burning will drive convec- 
tion, establishing an extended convective envelope with its 
base at slightly larger radii than the temperature peak. The 
object will likely grow to have a radius comparable to that 
of a giant star and correspondingly a relatively cool effec- 
tive temperature like the models presented in |Shen et al.| 
( |2012[ ). There are clear opportunities for future work in the 
self-consistent thermal evolution of these objects and their 
consequences for Type la supernovae, neutron stars, R Coro- 
nae Borealis stars and other phenomena. 
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APPENDIX A: VERIFICATION TESTS 

We perform a number of tests to confirm that our results 
are insensitive to the details of our approximations and nu- 
merical methods. A summary of these test runs is shown in 
Table |A1| Each run has an ID, which begins ZTn, where 
ZT represents "ZEUS testing" and n is an integer, indicat- 



ing that SPH simulation Pn of Dan et al. (20111 was used 



to generate the initial conditions. The results of these tests 
are discussed in the following sections. 

Al Resolution 

We confirm that our solutions are numerically converged by 
performing runs at different resolutions. We perform runs 
with 2/3 and 4/3 the linear resolution of the fiducial cal- 
culation. Fig. |A1| shows ID profiles from each of these runs 
after 10 4 s. There is only a small variation between the fidu- 
cial run (ZP5) and the high resolution run (ZT5-HR). The 
lower resolution run (ZT5-LR) also agrees quite well; the vis- 
ible variation is in the interior of the surviving WD, not the 
viscously evolving exterior. These results demonstrate that 
our simulations are converged in the quantities of interest. 
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ID 


Parameter 


Value 


ZT5-LR 


N r ,N e 


48,32 


ZT5-HR 


N r ,Ng 


96,64 


ZT5-alpha-m 


a 


icr 2 


ZT5-alpha-p 


a 


HT 1 


ZT5-hydro 


a-viscosity 


Off 


ZT5-visc-full 


a-viscosity 


All components 


ZT2-burn-apl3t 


Network 


aproxl3 


ZT2-burn-hecot 


Network 


HeCONe 


ZT2-burn-offt 


Network 


No Burning 


ZT5-IC1 


*end,SPH 


35 P 


ZT5-IC2 


*end,SPH 


35.6 P 


ZT5c-rnu-m 


Tv 


2.5 x 10 s cm 


ZT5c-rnu-p 




5.0 x 10 8 cm 



Table Al. Details of the test runs discussed in this appendix. 
ID is the run ID: ZT represents "ZEUS testing" and the num- 
ber indicates which initial conditions are being simulated. The 
string following the first dash briefly describes the parameter be- 
ing changed. Parameter is a description of the aspect of the run 
that was varied. Value is its value, 'this simulation had a lower 
resolution, N r ,Ng = 48,32 



A2 Independence of a 

We expect our simulations to be insensitive to the exact 
value of a so long as the hierarchy of time-scales in the 
problem remains unchanged. Specifically, a must not be so 
small that energy transport by radiation (or energy release 
from nuclear reactions) becomes important and not so large 
that the viscous time becomes less than the orbital time. 
Fig. |A2| shows that we observe only weak variation in our 
results with a in the range 0.01 — 0.1. The simulations are 
compared after a constant number of viscous times, such 
that atend = 3 x 10 2 s. 



A3 Viscosity Tensor 

Motivated by numerical simulations of the MRI we choose a 
prescription in which only the azimuthal components of the 
viscous shear tensor were retained. We relax this assump- 
tion and explore the effects of retaining all components of 
the tensor. This choice has virtually no effect on the evo- 
lution of the material near the temperature peak, as the 
large initial <f> velocity shear means that the non-azimuthal 
components of the tensor are small in comparison anyways. 
In the outer regions where azimuthal shear is not always 
so dominant, this choice can have an effect. In the fiducial 
case the evolution of the outer ~ Q.2Mq of material shows 
some minor differences. However, none of our conclusions 
are based on the detailed structure of the outer material, so 
this does not alter any of our conclusions. 



A4 Initial Conditions 

In order to confirm that our results are independent of small 
details of the initial conditions, we initialize our simulations 
with output from the same SPH calculation taken at three 
different times. By default, we start from the end of the SPH 
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Figure Al. The convergence of our simulations of the fiducial 
remnant with numerical resolution. The top panel shows ID tem- 
perature and entropy profiles and the bottom panel shows the 
ratio of the angular velocity to the Keplerian angular velocity. 
The overlap of the fiducial run (ZP5) and the high resolution run 
(ZT5-HR) indicate our simulations are converged in these quan- 
tities. 
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Figure A2. The variation of our simulations of the fiducial rem- 
nant with different values of a. The top panel shows ID tem- 
perature and entropy profiles and the bottom panel shows the 
ratio of the angular velocity to the Keplerian angular velocity. 
The profiles are shown after the same number of viscous times, 
at atend = 3 X 10 2 s. While there are small variations between 
runs, we see no significant change in our results for values of a 
spanning an order of magnitude. 
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A7 3D 

Moving to 3D makes the viscous evolution substantially 
more computationally demanding because of the strong 
timestep constraint imposed by our explicit evolution of the 
viscous terms, At v i ac ~ min((Ar) 2 /i/). A zone near the pole 
(6 « ir/(4Ne)) has size Ar = 2nr6/N 4) , where N$ is the 
number of (f> zones. The means that at the same r, 6 reso- 
lution, a 3D calculation will require evolving approximately 
N$ as many zones at timestep smaller by a factor of Ni. 
This issue can be helped somewhat by subcycling, that is 
advancing the viscous terms at A£ v i sc but integrating the 
rest of the hydrodynamics at A£cfl- 

In light of these issues, the simulation we perform is a 
simple one in which we initialize the same azimuthally sym- 
metric initial conditions used in 2D on a lower resolution 3D 
grid (N r = 48, Ng = 32, Afy = 32). We evolve the system for 
a substantially shorter time, only 1 x 10 2 s. Over that lim- 
ited time, we observe no qualitative differences which would 
affect our conclusions. 



Figure A3. The variation of our simulations of the 0.3 + I.IA/q 
remnant with different nuclear networks. The top panel shows ID 
temperature and entropy profiles and the bottom panel shows the 
ratio of the angular velocity to the Keplerian angular velocity. The 
results of our simple network (HeCONe) and the aproxl3 network 
are almost indistinguishable. We show the results of omitting the 
network entirely to illustrate the small impact of nuclear burning 
on the remnant structure over the viscous time-scales of interest. 



calculation, which for the fiducial model was after 35.7 ini- 
tial orbital periods had elapsed. (The secondary was tidally 
disrupted after 29 orbits.) The results we obtain with initial 
conditions from output taken 0.1 and 0.7 initial orbital peri- 
ods before the end of the calculation are virtually identical. 
The outcome of our calculation is insensitive to the dura- 
tion of the SPH simulation, so long as the remnant has had 
sufficient time to evolve towards axisymmetry. 



A5 Viscosity Cutoff 

As expected, we confirm that out results are insensitive to 
the location of the cutoff radius defined in Equation [6] and 
the surrounding discussion. 



A6 Nuclear Network 

As discussed in Section |2.3| we implement a simple 5 iso- 
tope (He,C,0,Ne,Mg) nuclear network. We confirm that this 
simple nuclear network reproduces the results of the more 
sophisticated aproxI3 network. Because of the high compu- 
tational cost of the full network, we perform these test calcu- 
lations at a lower resolution. We perform these tests on the 
0.3 + 1.1 He+CO system (ZP2) as it has the shortest burn- 
ing time of any of the He+CO mergers we consider. Fig. |A3| 
shows that the two networks give identical results. We also 
show the effect of omitting the nuclear burning, which does 
change the values of the thermodynamic quantities near the 
temperature peak, but does not alter the overall structure 
of the remnant. 
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